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Abstract 
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The equation v V u = studied by means of a compact 

finite difference scheme and numerical solutions are compared to the analytic 
inviscid (v = 0) solutions. The correct internal and external boundary 
layer behaviour is observed, due to an inherent feature of the scheme which 
automatically produces upwind differencing in inviscid regions and the correct 
viscous behaviour in viscous regions. 
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INTRODUCTION 


Consider a domain D with boundary T in the x y plane, in which 
u(v) is the solution of the time-independent convection diffusion equation 

f (u) + 8„(u) = V u (1.1) 

X y 

for u = u on F. Certain features of this problem can provide useful 

insights into the Navler-Stokes equations and their invlscid limit, the Euler 
equations (we employ terminology suggested by these applications in the 
following discussion) . In fluid dynamics it is common to obtain the time- 
independent solutions to these problems by solving the time dependent problem 
for large times. In this paper we Instead study (1.1) by a time-independent 
finite difference equation to which rather standard Iterative methods will be 
shown to apply. 

A large variety of finite difference methods can be employed to solve this 
problem in the sense that if u(v,h) Indicates the solution of a finite 

difference scheme in which h is a typical mesh length then u(v,h) u(v) 
as h 0 in, say, an L 2 norm. In particular, integrating (1.1) over a 
subdomain with boundary V' produces the conservative Integral 

formulation 

^ (f n + g n )ds = (^ V V u*n ds, (1*2) 

r- X y J,, 

T 

where n = (n ,n ) is the unit normal vector. Conservative finite 

— X y 

difference schemes approximate this equation in each cell. Nonconservative 
schemes instead approximate the equivalent differential equation 

2 

au + bu =vV u, 

X y ’ 


(1.3) 
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where a = 8f/3u and b = 8g/3u In each cell. Since u(v) is smooth both 
approaches produce the correct solution in the limit h ^ 0. 

The situation is more complicated for the "inviscid” weak solution u® 
which satisfies 

(f n + g n )ds = 0, (1.4) 

^ y 

for all closed curves T' . This is equivalent to the hyperbolic problem 


0 V 0 

a u + b u = 0, 
X y ’ 


(1.5) 


provided is differentiable. However, in the inviscid case boundary data 

may only be prescribed on part of P, say P^. Furthermore there may be 

nonunlcjue solutions with discontinuities (shocks) or closed characteristic 

2 2 

curves (recirculating flows) if a + b = 0 at points in D. 

Many of the essential facts concerning the relationship of the solution 
u(v) of (1-1) with the "physically relevant" solution u^ of (1.4) are 
suggested by formal singular perturbation arguments. Specifically, in the 
limit V ^ 0, we expect that the solution u(v) of (1.1) converges to the 
physical** weak solution u^ of (1.4) and boundary layers arise on the 
complement of P in P, i.e., that part of P on which no boundary data is 

required for the inviscid problem for u^. 

The use of finite difference schemes to calculate the solution u(v) in 
the limit v ^ 0 Introduces a second limiting process h ^ 0. It is natural 
to consider the different limiting paths by which u(v,h) ^ u^ as suggested 
by the following figure: 
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h 0 


u(v,h)- 


•u(V) 






V = V(h) 


V -»• 0 


Vh ->■ 0 




t 

u(h)— 


V 0 


h ^ 0 


Figure 1; Typical limiting processes considered In this paper. 


The solid arrows indicate limiting processes which are known to be 
convergent, assuming a consistent stable finite difference scheme. The dotted 
arrows denote limiting processes which in general may not be convergent, or 
may converge to a solution different from 

This paper is concerned with investigating a finite difference scheme 
which offers useful insights into this problem. The scheme has its origin in 
a compact exponential scheme described by Philips and Rose [1] for time— 
dependent problems; an application of this scheme to the Navier-Stokes 
equations was also described by Rose [3] . An effective time-independent 
solution method is provided by employing the flux-elimination technique 
described in Phillips and Rose [2]. The first few sections describe the 
application of these developments to equation (1.1). The relationship of the 
solution u(v,h) to u*^ is then explored by several numerical experiments 
and simple heuristic analyses. 



2. A COMPACT SCHEME 


Consider (!•!) in system form 


V + w 
X y 


+ 8y(u) 


W = V u 


In a square cell tt with sides of length 2h, , i, u, . i, , 

/2 i>J±V2 

'^i±V2>j* ^i,j±V2 average values of the analytic variables 

the corresponding sides of • 

J 



Figure 2: Location of finite difference variables in a cell 


7T 
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Because solutions of (2*1) are smooth, if h is sufficiently small any 
solution of (2.1) can be approximated locally by 

2 

vV u = au +bu 
X y 


where a,b indicate average values of f^, in the cell. This equation 

has elmentary solutions u = exp(ocx + 3y) where oi(a - va) + 3(b - v3) = 0. 
Any linear combination of four such solutions can be used to obtain a solution 
having the prescribed average values ^i±V2>j’ ^i>j±V2 sides of 

the corresponding average values Wj^^ are then related 

by four algebraic equations. The following algebraic system arises when a 
solution of the form 

U = Cl + C2 exp(^) + exp(^) + exp(^^- 


is considered. 


a) 6 v + 5 w = 6 f+5g 
^ X y X y ^ 


b) (Ujj - h p^ 6^)v = V 6^ u 


X X' 


c) fy -hp 6lw = v6 u 

'• y y’ y 


( 2 . 2 ) 


where 


d) h^(q^ V - 6^ w) = v[y^ - y^)u 


X X 


X y 


'x " ■ (“i+iA.j ■ 


\ “ ■ (“1+1/2.J 



etc*, and, if p = h/v, 0 = ap , 0 

X y 

= bp , then 


p(0) = coth 0 - 0 ^ 


111 

q(0) = 0"^ p(0) 

"x ^ "('x’’ 


The simple behaviour of p(0), q(0) 

are shown in Figure 3 


Figure 3. Graph of p(0), q(0). 


The dimensionless parameters 0 >0 are the cell Reynolds numbers and 

X y 

provide a measure of the relative importance of diffusion versus transport 
effects in a cell* An important feature of this scheme is the manner in which 
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upwind differencing arises when l0| °o. From (2.3), (2.4), 


p(6) -- 

sgn(6) 

lei 

large 


e/3 

lei 

small 

q(9) 

ler^ 

lei 

large 


1/3 

lei 

small. 


In (2.2b) (p^ - h represents a weighted average of the values 

v^^ j the two cell sides. For |0 | -v «> the weights tend to 0 or 1 

thereby selecting one or another of the values depending on the sign of 0. 
For l0| 0 the weights tend to V 2 giving a simple central average v. 

An important point to note is that equation (2.2a) can be interpreted as a 
finite difference form of the conservative integral equation (1.2). Summation 
over all the cells in domain D yields 

^ (f n + g n )h = ^ (v n + w n )h, (2.6) 

p ^ y p X y 

so the scheme is globally conservative. This is Important in ensuring the 

correct calculation of nonlinear shocks, as will be shown later. 

The scheme (2.2) is called compact because it involves only the values 
of u, V, w associated with an Individual cell. Any relationships between 
values in other cells arise from imposing the further condition that values 
are continuous across cell boundaries. When a and b are constant in all 
cells an energy estimate similar to that in [1], [2] can be obtained from 

(2.2); the result is that the solution u(v,h) converges as h 0 to a 

p 

smooth solution u(v) of (2.1) with O(h^) accuracy. 
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3. FLUX ELIMINATION 
Let 


V = [u v,A v,y w,A w] 
X * X y y ^ 


U = fy u, A u, y u, A u) 
^ X ’ X ’ y " y ^ 


(3.1) 


where A = h 6. Equations (2.2) may be solved for V to obtain 


V = R(U) 


(3.2) 


where 


R(U) = a 


-1 


oA u + p[(y -y)u + pq(A f+A e)l 
X ^x*- X y ^y X y ^ ^ 


(y - y )u + p q (A f+A g) 
X y ^y X y ^ 


oA u+p[(y -y)u+pq(A f+A g)l 

y y*- y x ^x^ x y 


(3.3) 


(y - y )u + p q (A f+A g) 
y x^ ^ ^x^ X y 


in which a = q + q . 

Following the method outlined in Phillips and Rose [2] consider two 
neighboring cells and j having the common values • and 

Vi+ I /2 , j associated with their common side. Clearly 


"i+Vz-J ■ ("x ■ (■': 




SO that 
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(1 1 0 O)r^^j(U) = (1 -1 0 0)^1+!, j(U). 

Similarly considering the value i+ V 2 to cells , and 

^i j+1 also have 

(0 0 1 1)^1, j(U) = (0 0 1 j+1 (U). 


Using (3.3) there results 

a) P ^ 

“ ^X " Py^*^ P ^ + ^y 8)]} =0 

(3.4) 

b) yy{o"^[(My - Mx^" ■*■ P ^^X^^x ^ ■'■ ^y 

“ \^^y " + Py P ^[(Py " Px^'^ ■*■ P ‘’x^^x ^ ■'■ ^y s)]} = 0- 

The values of u related by each of these equations are indicated in 
Figure 4. We refer to (3.4) as the flux-eliminated form of (2.2). From its 
solution V and w may be calculated from (3.2). 

When only u is prescribed on the boundary then (3.4) applies as 
indicated at all interior points. When, more generally, v or w are 
prescribed at boundary points additional equations expressing these prescribed 
values are obtained from (3.2). Finally we call attention to the fact that 
the parameters h and v occur in (3.4) only in the combination p = h/v. 
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(a) (b) 

Figure 4: Stencils for (3.4a,b) 

4. LIMITING FORMS OF FLUX-ELIMINATED EQUATIONS 

As stated in the Introduction we shall be Interested in studying the 
solution u(v,h) of (3.4) as h 0, v(h) 0. Depending on the function 

v(h), as h 0, p = h/v approaches 0, a constant, or In this section 

we shall describe some results of formally applying the limits p -»■ 0,“ in 
(3.4) using the asymptotic values for p,q given in (2.5). For our purpose 
it will be sufficient to examine only (3.4a). 

Case 1 : p ^ 0 

When 6 2 0 then p 2 *12 (3.4a) reduces to 

*“i- V2 .J “1+V2 f] “1+3/2, j' 

• ^(“l.j-1/2 “l.j+1/2 * “1+1.J-V2 + “l+l.J+Vz) - “• 
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Case 2 ; p -»■ » 

When |6| » 1 then p sgn(6), q ~ l9| ^ and so in this case different 
limiting forms arise depending on sgn(a), sgn(b) where a = f^, b = g^. We 
write 

" ®i.j " ^i+l,j 

" ^i,j " Vl.j 

(i) b^>0 

Equation (3.4a) reduces to 

fb - a lu. 1 , . + (b + a lu,, 2b u, . i, '^0. 

The relevant stencil is shown below with arrows indicating the direction of 
the inviscid characteristic. Note that the reduced equation uses only 
"upwind” information to calculate u^^l ^2 j* particular if these two cells 

are the entire computational domain then the prescribed boundary values on the 
upwind side of u^^ 1^^ j solely determine its values. This corresponds to 
the analytic situation in which the inviscid weak solution requires boundary 
data only on the 'inflow' part of the boundary at which characteristics enter 
the domain. 



Figure 5: case (i) 
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(ii) < 0, a_ > 0, > 0 (Ibl » lal) 

Equation (3.4a) reduces to 


" 1 - 1/2 ,j 1 <»- - =+>" 1 + 1/2 ,j - =+ " 1 + 3 / 2 , j 


-2a u. . 1 / + 2a, u. ,, . i. ~ 

- i.J-V2 + i+l,j-V2 - 


~ 0 . 


Again u^^l ^2 is influenced only by upwind points as shown below. 



Figure 6: case (il) 

(iii) a^ > 0, a_ < 0, b^ > 0 (Ibl » lal) 

Here 


'^i+V2,j 2 ^'^i,j-V2 '^i+l , j- V 2 ^ - 


0 . 



Figure 7: case (iii) 
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In interpreting these observations it must be remembered that 
corresponding equations arise from (3.4b). These examples provide useful 
insights in understanding the test problems considered in the next section. 
In these problems p varies from 0.3 to 40 with corresponding values of 6 
ranging from 0 to 40 so the actual finite difference equations can be very 
close to the different limiting forms in different parts of the domain. 


5. NUMERICAL EXAMPLES 

We have indicated earlier that we may expect that the compact scheme (2.2) 
yields a solution u(v,h) which converges to the solution u(v) of (2.1) 

o 

with an L 2 error of order h (c.f. Example 1, below). In this section we 
are primarily interested in studying the L 2 error E between u(v,h) and 
the inviscid solution u^ as h 0, v 0, by means of numerical 
experiments. 

The first question which arises concerns the possibility of solving (3.4) 
for u(v,h) when p ->- 00 . Although several Iterative methods were studied 
with, generally, successful results our report will concentrate on the use of 
the Gauss-Seidel method. One reason for our doing so stems from favorable 
results, as yet unpublished, by our colleague T. N. Phillips concerning the 
treatment of systems closely related to (3.4) by multigrid methods using 
Gauss-Seidel as the underlying iterative scheme. 

A central question motivating this study is how effectively internal and 
external boundary layers arising from approximating u^ by u(v,h) can be 
localized. A related question concerns the relative influence of the 
viscosity v compared to Inaccuracies arising from the fact that the scheme 
has truncation errors of order h^. In order to help discuss this we 
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2 

introduce a parameter x = h /v which together with p = h/v, will be used 

to study the solution as h -»• 0, v(h) -»■ 0. 

Another issue concerns new insights into upwind difference schemes. As 
discussed earlier |6l « 1 leads to a diffusion-type limit of (3.4) while 
|0| » 1 leads to an inviscid upwind-type treatment. Both can occur if 
P 1 and a (or b) passes through zero in some region. The functions 
p,q in (3.4) automatically handle this transition with a "viscous” treatment 
near a 0 and an "almost-inviscid" treatment elsewhere. 

The final matter concerns the question; under what limiting processes do 
nonconservative finite difference schemes yield the correct as 3 nnptotic 
inviscid result when h > 0, v(h) -*■ 0? 

The examples which follow address these questions. The first two examples 
graphically Illustrate the behaviour of u(v,h) in cases in which the 
inviscid limit u® is either linear or nonlinear. In the next five examples 

the L 2 error E between the solution u(v,h) and the inviscid, analytic 

solution u® Is computed for different values of v,h and the results are 
compared using simple heuristic analyses. The final example calculates the 
1^2 error for a problem using a nonconservative form of (3.4). 

Example 1 

f and g are linear and are defined by 

f = au a = x-V 2 

g = bu b = l- y. 

Note that f^ + g^ = a u^ + b u^ since a^ + b^ = 0. Hence the inviscid 
solution u® is constant along characteristics defined parametrically by 
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dx 

ds 


a, 



i.e., (1 - y)(V 2 “ x) = constant. 


The domain considered is 0 < x < 1, 0 < y < 1 and so the only inflow 
boundary on which boundary data needs to be specified is on y ~ 0. The 
boundary values are chosen to be 


u^(x,0) 


sin Ti ( V 2 - x) , 


so that the solution on the interior is 

u^(x,y) = sin[TT(l - y)(V2“ x)] . 


The boundary conditions for the finite difference problem are chosen to be 

u(x,l) = u^(x,l), u(0,y) = -1, u(l,y) = 1 

SO that boundary layers arise at x = 0,1. 

Figure 8 shows contour plots of u(v,h) for several different values of 
V and h. (a), (b) , (c) illustrate ordinary convergence arising from a 

sequence in which v is kept fixed and h is reduced by factor V 2 each 

time; little visual difference appears. (a), (d), (e) are a sequence in 

which p = h/v is kept fixed and h is reduced by factor 2 each time; 

observe that the width of the boundary layers decreases. Finally, (a), (f), 
(g) illustrate, a sequence in which h is kept fixed and v is decreased by 
a factor 4 each time. Note that the boundary layers at x = 0,1 become 
smaller until they are confined to one interior cell. 



(a) 

V = 1/5, N = 10 


(d) 

V = 1/10, N = 20 


(e) 

V = 1/20, N = 40 



V = 1/20, N = 10 V = 1/80, N = 10 

Figure 8. Contour plots for Example 1; (a), (b), (c) illustrate 
ordinary convergence (h -»• 0, v = const.); in (a), (d), 

(e) p = const., h ->• 0; in (a), (f), (g), h = const., 

V ^ 0. (In (f), (g) boundary layer thickness is confined 
to one cell.) 
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Example 2 

This nonlinear example was suggested by our colleague E. Tadmor to test 
the ability of the scheme to converge to the correct "physical" inviscid 
solution for a problem with an infinite range of formal inviscid solutions. 
We consider 

- 2 3 

f = u , g = u . 


The test case has an expansion fan attached to a shock 


u^(x,y) 


0.8 


2y 

3(x - 0.3) 


X 


-0.3 

y 


< 


1 

1.2 


1 X - 0.3 1 

1.2 y 0.6 • 


V -0.2 


1 X - 0.3 
0.6 y 


The boundary conditions for the finite difference solution are 
u(v,h) = u^ on all four sides so that the only sharp gradients occur at the 
shock. 

Figure 9 shows contour plots of u(v,h) for fixed P and decreasing h. 
In (a) the shock is indistinguishable from the expansion fan but in (c) the 
difference is obvious, and it is also clear that the numerical solution is 
converging to the physical inviscid solution. 
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(a) 


(b) 

2 ^ 

Figure 9» Example 2, f = u , g = u*^; p 


(c) 


const. , h 0. 


Example 3 

This linear example uses the same definition of f and g as Example 1, 
and has the same inviscid solution u^. However this time the finite 

difference scheme has boundary conditions u = u^ on all four sides of the 
square domain. Three cases were run, one with a 20 x 20 grid, and two with 
80 X 80 grids keeping p = h/v fixed in one case, and x = h /v fixed in the 
other. The table below shows the (i) number of Gauss-Seidel iterations, 
(starting from initial conditions u = 0) ; (ii) the L 2 error E; (iii) E 
divided by the error E* for the 20 x 20 grid, and (iv) the predicted ratio 
E/E based on the following simple analysis. 

Let Lq be the inviscid differential operator and L(v,h) be the viscous 
finite difference operator. Then 

Lq u^ = 0 and L(v,h) u(v,h) = 0, 
so 

2 

L('^»h)(u - u(v,h)) = (L(v,h) - ~ h^ • function (~) 
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and 


0 


u 


u(v,h) =0 on boundary; thus 


llu° - u(v,h)ll^ ~ e(T) as h > 0, 

^2 

where e(T) is some function of the truncation error due to viscosity* 


Table !• Results for Example 3 



N = 20 
P = 10 
T = 0.5 

N = 

p = 10 

T = 0.125 

80 

P = 40 
T = 0.5 

No. of Iterations 

15 

55 

25 

E 

.334 (-2) 

.821 (-3) 

.216 (-3) 

E/E* 

1 

0.245 

0.0647 

Predicted Ratio 

1 

0.25 

0.0625 


Example 4 

This example has linear f,g 

f = 0.5 u, g = u 

with discontinuous boundary data for which the inviscid solution u^ is 

Q ( 0.8 X < 0.3 + 0.5y 

u = 

(-0.4 X > 0.3 + 0.5y 


The analysis for this case is not easy but leads to the result 



where e is a function with asymptotic behaviour. 


T > 0 



===> 


E 



Thus when t is small the error E is due not to the second order truncation 
but due rather to the viscous smearing of the contact discontinuity. 


Table 

11 • Results 

for Example 4 



N = 20 

N = 

80 


P = 10 

P = 10 

P = 40 


T = 0.5 

X = 0.125 

X = 0.5 

No. of Iterations 

10 

15 

15 

E 

.144 

.0969 

.0901 

E/E* 

1 

0.673 

0.625 

Predicted Ratio 

1 

0.707 

0.5 


Example 5 

This example is a two-dimensional version of Berger's equation with a 
shock: 

f = 0.5 u^, g = 0.4 u 


and the inviscid solution 


u 


is 


0.8 


21 


0 

u 


-0.4 


X < 0.3 + 0.5y 
X > 0.3 + 0.5y 


In the finite difference solution the shock becomes an internal boundary 
layer with width of order v. If p = h/v is kept fixed as h 0, this 
shock layer is spread over the same number of cells and hence 


E ~ h 




e(p). 


As p 0 , e ~ p~ ^^2 ===> E ~ V ^^2 gg when P is small the error is due 

solely to a physically well-resolved shock. 


Table 

III. Results 

for Example 5 



N = 20 

N = 

80 


p = 10 

P = 10 

p = 40 


T = 0.5 

X = 0.125 

X = 0.5 

No. of Iterations 

35 

85 

85 

E 

.102 

.050 

.048 

E/E* 

1 

.49 

.47 

Predicted Ratio 

1 

.5 

.5 


Example 6 

This example is exactly the same as Example 1 except for different values 
of h and v. The analysis cannot be expressed in one simple equation. If 
p ^ 0 and h 0 then E v ^2 since the boundary layer has thickness of 
order v and so it is well-resolved. If p is kept fixed the boundary layer 
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spans a fixed number of cells and so E ~ h '^2 , j j^gp^ fixed it has 

an 0(h) effect on just the interior points adjacent to the boundary and so 
^ u3/2 


Table IV- Results for Example 6 



N = 20 
P = 10 
T = 0.5 

N = 

P = 10 
T = 0.125 

80 

P = 40 

T = 0.5 

No- of Iterations 

15 

55 

25 

E 

.448 (-1) 

.254 (-1) 

.51 (-2) 

E/E* 

1 

0.567 

0.114 

Predicted Ratio 

1 

0.5 

0.125 


Example 7 

This example is the same of Example 3 except that the boundary condition 
at y = 1 is u(x,l) = 1 which produces a boundary layer at y = 1- This 
boundary layer is different from the ones in Example 6. In Example 6 the 
boundaries were 'outflow^ boundaries, i-e«, the inviscid characteristics were 
pointing outwards across the boundary- In this example the characteristics 
are tangential to the boundary at y = 1, similar to a stagnation point flow 
in fluid dynamics • 

The physical boundary has thickness of order v ^2 g^ j ^g j^gp^ 

fixed as h 0 it spans a fixed number of cells and E ~ h ^^2 . Hence 


E ~ h ^'*2g(T) 
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and as x 0, e 


x“ ===> E ~ . 


Table V. Results for Example 7 



N = 20 
P = 10 
X = 0.5 

N = 

P = 10 
X = 0.125 

80 

P = 40 
X = 0.5 

No. of Iterations 

20 

60 

25 

E 

.149 

.107 

.074 

E/E* 

1 

0.718 

0.498 

Predicted Ratio 

1 

0.707 

0.5 


Table VI summarizes the results of Examples 3-7. 

Table VI. Summary of Asymptotic Behaviour of E 


E 

V = const p = const 

X = const 

Smooth u® 

V h 

h2 

Linear, discontinuous u® 

V h ^^^4 

h^/2 

Nonlinear, discontinuous u® 

V h 

h^/2 

Smooth u^, “wrong** outflow 

boundary conditions 

V /2 h^/2 

^3/2 

Smooth u^, **wrong** tangent 

boundary conditions 

V h 

hV2 
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Example 8 

This example studies a nonconservative form of the compact equations 

obtained by replacing the term 6 f + 6 g in (2.2a) by a6u + b6u 

X y X y 

where a, b are average values of f^ and g^ in the cell. The test case 
is 

f = u^, g = 0.4 u 


with an inviscid solution 

Q I 0.8 X < 0.3 + 0.5y 

u = 

1-0.4 X > 0.3 + 0.5y 

The error analysis shows that E has two components, E^ due to shock 
smearing, and E 2 due to an incorrect shock angle. 

E^ (shock width) 

~ v^/2 = h^^2p-V2 


^2 ^ (conservation loss per unit shock length) 


V2 


[ (// cells across shock) x (// cells along shock) 
(truncation error per cell)] 

^ [ V X i X lii] 


~ P 


V2 


25 


The two error components have the same order of magnitude when p 

^2 T7 nof*AC an A If n ^ ^2 ^ 


h ^/2 


If p « h 


, dominates and. If p » 


dominates. Hence, 


E ~ 


M 

V2 


p « h 
p » h 


V2 

V2 


V9 

and at each fixed h, E is minimized by choosing p h ^ , and 

E ^ h . 
min 

Table Vila shows results from several experiments and Table Vllb has the 

V4 . /i, V2 


corrresponding values of ^min^^ 


and 


which, according to the 


analysis, should be constants. 


Table Vila: 


E 

P 

20 

N 

80 

0.31 

.267 


.162 

0.62 

.220 


.121 

1.25 

.168 


.133 

2.5 

.228 


.266 

5 

.376 


.393 


Table Vllb 


N 



20 

80 


.355 

.362 


.18 

• 

00 
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6. CONCLUDING REMARKS 

The experiments reported in the previous section suggest that the compact 
scheme (2*2) provides an effective means of approximating both u(v) and u® 
and that the relationship between these solutions which are suggested by 
singular perturbation arguments is maintained, as indicated by Figure 1, by 
the finite difference scheme as well* A fact of potential practical 
importance is that boundary layers can be confined to a single computational 
cell* This feature is a relevant factor as well in selecting the conservative 
form of (2*2) instead of the nonconservative form (we remark in passing that 
it would have been preferable to have treated the terms p6, q6 in (2*2) by 
6p, 6q instead)* 

Finally, the fact that the flux-eliminated equations (3*4) can be treated 
by a Gauss-Seidel iterative method indicates that a wide variety of more rapid 
iterative methods can also be employed* In contrast, the use of time-stepping 
methods to obtain u(v) or u® can be seen to result in a more slowly 
convergent Jacobi-type iteration scheme* 
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